home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / RANDOMTP / RAND.PAS next >
Pascal/Delphi Source File  |  1990-01-15  |  7KB  |  192 lines
  1. {$A+,B-,D-,E-,F-,I-,L-,N-,O-,R-,S-,V+}
  2. Unit Rand;
  3. Interface
  4. { RAND.PAS - This unit implements the "minimal standard" random number
  5.   generator as described in the October '88 "Communications of the acm"
  6.   article by Park and Miller.  The algorithm is coded in Inline assembler
  7.   using a variation of the technique described in the January, 1990
  8.   CACM article by Carta.  The floating point result returned by function
  9.   Random differs from the Park and Miller function only in the interval
  10.   bounds: [0.0 <= Random < 1.0) for this unit as opposed to
  11.   (0.0 < Random < 1.0) for Park/Miller.  The IRand Inline macro (see below)
  12.   has been verified to yield the same seed value as the Park/Miller function
  13.   after 3 million iterations, including several "overflow" conditions
  14.   (see Carta).
  15.  
  16.   To eliminate the "parallel [hyper]planes" pattern in supercritical
  17.   applications, a variation is provided that "Shuffles" the output of the
  18.   minimal standard generator as described in Knuth, "The Art of Computer
  19.   Programming", Vol 2 (Seminumerical Algorithms), page 32, Algorithm "B".
  20.  
  21.   Also provided is a routine to assist in selection of an unbiased sample
  22.   from a population of known size.
  23.  
  24.   No attempt was made to do a statistical analysis of the output of these
  25.   generators.  However, their performance has been analyzed extensively in
  26.   the works cited above and their references.
  27.  
  28.   This work is contributed to the public domain by the author:
  29.  
  30.   Edwin T. Floyd
  31.   #9 Adams Park Court
  32.   Columbus, GA 31909
  33.   (404) 322-0076 (home)
  34.   (404) 576-3305 (work)
  35.  
  36.   First release, 1-15-90.
  37. }
  38. Const
  39.   RandSeed : LongInt = 1;{ Random number seed: 0 < RandSeed < RandModulo }
  40.   Shuffled : Boolean = False; { True if shuffle table is full }
  41.   ShuffleBits = 5;       { 5 => 32-way shuffle, 6 => 64-way, etc }
  42.   RandMult = 16807;      { =7**5; RandMult must be expressable in 16 bits.
  43.                            48271 may give better "randomness" (see ACM ref.),
  44.                            but shuffling is much more effective (see Knuth). }
  45.   RandModulo = $7FFFFFFF;{ =2**31-1 = 2147483647 }
  46.  
  47. Function IRand : LongInt;
  48. { Return next "minimal standard", 31 bit pseudo-random integer.  This function
  49.   actually computes (RandSeed * RandMult) Mod (2**31-1) where RandMult is
  50.   a 16 bit quantity and RandSeed is 32 bits (See Carta, CACM 1/90).  Apart
  51.   from the load/store of RandSeed, this may be the fastest possible
  52.   implementation of this function with standard 8086 instructions on 16 bit
  53.   registers. (The glove is down!) :-}
  54. Inline(
  55.   $A1/>RandSeed+2/       {         mov     ax,[>RandSeed+2]}
  56.   $BF/>RandMult/         {         mov     di,>RandMult}
  57.   $F7/$E7/               {         mul     di}
  58.   $89/$C3/               {         mov     bx,ax}
  59.   $89/$D1/               {         mov     cx,dx}
  60.   $A1/>RandSeed/         {         mov     ax,[>RandSeed]}
  61.   $F7/$E7/               {         mul     di}
  62.   $01/$DA/               {         add     dx,bx}
  63.   $83/$D1/$00/           {         adc     cx,0 ; cx:dx:ax = Seed * Mult }
  64.   $D0/$E6/               {         shl     dh,1 ; split p & q at 31 bits }
  65.   $D1/$D1/               {         rcl     cx,1}
  66.   $D0/$EE/               {         shr     dh,1 ; cx = p, dx:ax = q }
  67.   $01/$C8/               {         add     ax,cx}
  68.   $83/$D2/$00/           {         adc     dx,0 ; dx:ax = p + q }
  69.   $71/$09/               {         jno     done}
  70.   $05/$01/$00/           {         add     ax,1 ; overflow, inc(p + q) }
  71.   $83/$D2/$00/           {         adc     dx,0}
  72.   $80/$E6/$7F/           {         and     dh,$7F ; limit to 31 bits }
  73.                          {done:}
  74.   $A3/>RandSeed/         {         mov     [>RandSeed],ax}
  75.   $89/$16/>RandSeed+2);  {         mov     [>RandSeed+2],dx}
  76.  
  77. Procedure Randomize;
  78. { Initialize RandSeed using BIOS clock tick counter }
  79.  
  80. Function Random : Real;
  81. { Return pseudo-random, floating point number uniformly distributed on
  82.   the interval [0.0 <= Random < 1.0) }
  83.  
  84. Function SIRand : LongInt;
  85. { Return shuffled, 31 bit, pseudo-random integer }
  86.  
  87. Function SRandom : Real;
  88. { Return shuffled, pseudo-random, floating point number uniformly distributed
  89.   on the interval [0.0 <= SRandom < 1.0) }
  90.  
  91. Function SelectRecord(Sample, Population : LongInt;
  92.   Var Selected, Index : LongInt) : Boolean;
  93. { Use this function to select an unbiased, random sample from a population of
  94.   known size, without replacement. Set Index := 0 before the first reference to
  95.   SelectRecord and test each record in the population for membership in
  96.   sample via reference to SelectRecord with constant Sample <= Population
  97.   parameters.  SelectRecord is guaranteed to return True for exactly "Sample"
  98.   calls, randomly distributed over the population.  This routine will re-select
  99.   the same sample given the same initial RandSeed and shuffle table state.
  100.   For a discussion of the algorithm, see Knuth, "Art of Computer Programming",
  101.   Vol 2 (Seminumerical Algorithms), page 137. Contrived example:
  102.  
  103.   Program DealOnePokerHand;
  104.   Uses Rand;
  105.   Var
  106.     i, j, k : LongInt;
  107.     CardFile : Text;
  108.     Card : String;
  109.  
  110.   Begin
  111.     Assign(CardFile, 'CARDNAME.TXT');
  112.     Reset(CardFile);
  113.     k := 0; (Start with zero index)
  114.     Randomize;
  115.     For i := 1 to 52 Do Begin
  116.       ReadLn(CardFile, Card);
  117.       If SelectRecord(5, 52, j, k) Then WriteLn(Card);
  118.       (Select unbiased sample of 5 from a population of 52)
  119.     End;
  120.   End.
  121. }
  122.  
  123. Implementation
  124. Const
  125.   RandomizeCallCount : Word = 1;
  126.   ShuffleShift = 16 - ShuffleBits;
  127.   ShufTableEnd = $FFFF Shr ShuffleShift;
  128.  
  129. Var
  130.   ShufTable : Array[0..ShufTableEnd] Of LongInt;
  131.   NextOut : Word;
  132.  
  133. Procedure Randomize;
  134. Var
  135.   BiosClock : LongInt Absolute $40:$6C;
  136.   ClockLow : Byte Absolute $40:$6C;
  137.   i : Byte;
  138. Begin
  139.   RandSeed := BiosClock And RandModulo;
  140.   If RandSeed < 1 Then RandSeed := 1;
  141.   If RandSeed >= RandModulo Then RandSeed := Pred(RandModulo);
  142.   Shuffled := False;
  143.   Inc(RandomizeCallCount);
  144.   For i := 1 To (RandomizeCallCount + ClockLow) And $FF Do
  145.     If Boolean(IRand) Then;
  146. End;
  147.  
  148. Function Random : Real;
  149. Begin
  150.   Random := Pred(IRand) / Pred(RandModulo);
  151. End;
  152.  
  153. Function SIRand : LongInt;
  154. Var
  155.   y : LongInt;
  156. Begin
  157.   If Not Shuffled Then Begin
  158.     For NextOut := 0 To ShufTableEnd Do ShufTable[NextOut] := IRand;
  159.     y := IRand;
  160.     NextOut := Word(y) Shr ShuffleShift;
  161.     Shuffled := True;
  162.   End;
  163.   y := ShufTable[NextOut];
  164.   ShufTable[NextOut] := IRand;
  165.   NextOut := Word(y) Shr ShuffleShift;
  166.   SIRand := y;
  167. End;
  168.  
  169. Function SRandom : Real;
  170. Begin
  171.   SRandom := Pred(SIRand) / Pred(RandModulo);
  172. End;
  173.  
  174. Function SelectRecord(Sample, Population : LongInt;
  175.   Var Selected, Index : LongInt) : Boolean;
  176. Begin
  177.   If Index = 0 Then Begin
  178.     Selected := 0;
  179.     If Sample > Population Then Begin
  180.       WriteLn('Sample must be <= Population in RAND.SelectRecord');
  181.       Halt(1);
  182.     End;
  183.   End;
  184.   If (SRandom * (Population - Index) < (Sample - Selected))
  185.   And (Selected < Sample) Then Begin
  186.     SelectRecord := True;
  187.     Inc(Selected);
  188.   End Else SelectRecord := False;
  189.   Inc(Index);
  190. End;
  191.  
  192. End.